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Abstract 

We consider a class of (2+1) dimensional baby Skyrme models with potentials that have more 
than one vacum. These potentials are generalisations of old and new baby Skyrme models; they 
involve more complicated dependence on (j)^. We find that when the potential is invariant under 
03 — > —03 the configurations corresponding to the baby skyrmions lying "on top of each other" 
are the minima of the energy. However, when the potential breaks this symmetry the lowest field 
configurations correspond to separated baby skyrmions. We compute the energy distributions for 
skyrmions of degrees between one and eight and discuss their geometrical shapes and binding 
energies. We also compare the 2-skyrmion states for these potentials. Most of our work has 
been performed numerically with the model being formulated in terms of three real scalar fields 
(satisfying one constraint). 



1 Introduction 

The Skyrme model describes a non-linear theory for SU(2) valued fields which has soliton solutions. 
Although the potential term is optional in the (3+1) dimensional nuclear Skyrme model its presence 
is necessary in (2+1) dimensions to ensure the stability of these solitonic solutions. In this paper 
we discuss multisolitons in a two-dimensional version of the model for a class of potentials which 
generalises the previously studied cases. In the literature we can find three specific potentials that 
have been studied in some detail. 

The holomorphic model has only one stable solution (describing one skyrmion) and this solution 
has a simple analytical form jj] . The old baby Skyrme model has stable solutions for any number 
of skyrmions. In they were shown to lead to a crystalline lattice of skyrmions. The new baby 
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Skyrme model was studied in detail in |^ where it was shown that its solutions corresponded to field 
configurations with radially symmetrical energy densities - which correspond to many skyrmions 
lying "on top of each other" . We generalise these studies by considering potentials with a more 
general dependence on 03 thus leading to models with more vacua. 



2 A Skyrme model for a class of potentials 

The Skyrme model was first proposed by T.H.R Skyrme 0| in 1960. Its classical solutions fall 
into various classes characterised by a topological number to be identified with the baryon number. 
Thus the model can describe mesons as well as various baryon configurations. The baby Skyrme 
models [Q are (2+1) dimensional versions of the Skyrme model and its topologically nontrivial 
solutions are called baby skyrmions. 

The Lagrangian density of the model contains three terms, from left to right: the pure S'^ sigma 
model, the Skyrme and the potential terms: 



The vector (p is restricted to lie on a unit sphere hence (p ■ <j) ~ 1 . 



Note that to have finite potential energy the field at spatial infinity cannot depend on the polar 
angle 9. 

lim $ir,9) = (2) 

r — *oo 

Hence this boundary condition defines an one-point compatification of i?2, allowing us to con- 
sider (p on the extended plane -R2IJ00 topologically equivalent to S^. In consequence, the field 
configurations are maps 

M:S^ — >S^. (3) 
which can be labeled by an integer valued topological index Q: 

Q = ^e"" / dxdycj^a {d^cPb) {dy^c) . (4) 

As a result of this non-trivial mapping the model has topologically nontrivial solutions which 
describe "extended structures" , namely, baby skyrmions . The different choices of the potential 
term lead to various shapes of the energy density of these baby skyrmions. 

The equation of motion for a general potential that depends on ^3 takes the form: 



~{dj-dj){d''$- d^$)4>a] + lSv^{5a3 - <Pa4>3) = 0. (5) 

It is convenient to rewrite this equation as: 

du(pa = K^bTb {$, dt$, diipj (6) 

with 

Kab - (1 + 29sd,$ ■ d,$)5ab ~ 29sd,(t>ad,^b (7) 
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and a rather complicated expression for Tb- 

Then to find a solution to this equation we invert the matrix K and simulate the time evolution 
by a 4th order Runge Kutta method supplemented by the imposition of a correction due to the 
constraint (p ■ cj) ~ 1. 

The kinetic and potential energy densities of the baby Skyrme model are: 

]C^{dt$-dt$)(l + 29 s (d,$- dS) ) - 26 s (dt$ ■ d^$) ' (8) 

V = (d,$- dS) + Os • ' - (a,0 • 9,0) (a^ • d,$) + OvVih). (9) 

The potentials we want to consider in this paper are of the form 

V{h) = (1 - 03)G(03), (10) 

where G{cj)3) = 1 or (1 + ^3) or 0| or (1 + (j>3)(l)^. Thus this potential vanishes at (^3 = 1 and also 
at 4>3 = —1 and/or ^3 = 0. 

The total energy takes the form: 

E = J dxdy (V + /C) . 

3 Static Solutions 

In this paper we are primarily interested in finding minimal energy configurations corresponding 
to many skyrmions. We want to see, for different forms of the potential, what such configurations 
are like and what properties they have. 

We will consider a numerically found configuration to be a multiskyrmion if it is satisfies certain 
numerical checks for a local minimum and if, moreover, its energy satisfies 

En < Ek + El for all integers 1 < £ and k < n such that k + £ = n 

We have included this condition in our definition because we want multisoliton to be stable with 
respect to decay into multiskyrmions of smaller degree. 

There are many ways of finding static solutions. First, as the potential is only a function of 
(f)3, there is a symmetry corresponding to rotations around the (jjs axis. So choosing the spatial 
dependence conveniently it is clear that fields which correspond to a generalised multihedgehog 
ansatz will be static solutions of the equations of motion . Such fields are given by ie 

' sin[/(r)] cos{n9 — x) \ 

sin[/(r)] sm{ne - x) ■ (11) 
cos[/(r)] / 

where (r,0) are polar coordinates in the a; — y-plane and the function /(r), called the profile function, 
is required to satisfy certain boundary conditions to be specified below. The angle x is arbitrary, 
but fields with different x are related by an iso-rotation and are therefore degenerate in energy, n 
is a non-zero integer and equals to the topological charge. 

For the field be regular at the origin the profile function f{r) has to be satisfy : 

/(O) = mn, (12) 

where m is an integer. 
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We choose the vacuum at infinity to be (fP — (0, 0, 1) and this results in another boundary 
condition, namely: 

lim J{r) = 0. (13) 

r — >oo 

The total energy of the field configuration then takes the form: 
1 p /^y2 2sinV, 



^"^^^^2io ""^''V +n'—j^{l + 2esf) + 9vY{f)j, (14) 

where f — ^ and V(/) = V{(f>3). To determine the profile function f(r) we treat ( p^ as a 
variational problem and we get a second-order ODE for /: 

^ 29sn^sufj\ f" ^Osn^sin^f _^ 20 sn^ sin f cos ff' \ ^, 

r J \ r"^ r J 

n2 sin /cos/ 9v dN{f) 
r '^T^ = ^- 

This equation then has to be solved numerically (ie via the shooting method). 
The topological charge takes the form: 

Q^-'^j^drl^-^^ -^[cos/(oo)-cos/(0)]. (16) 

This equation shows that to have an integer value for the topological charge, m in (^2|) must 
be an odd number. In this paper we consider the solutions with m = 1. 

The behaviour of solutions of this equation near the origin and for large r can be deduced 
analytically and was also discussed in The result is that 

• For small r, the profile function behaves as 

/ ~ TT + C„r" (17) 

and so 



/ ~nC„r"-i (18) 

df 



as long as '^^}P tends to zero at this point 



At large r, the ODE reduces to 



1 +-J-:^-l-^ ^llarger-"- (.1^) 

Using the boundary conditions, we see that if y^ljarge r ^ "'^ ^^^^ profile function decreases 
exponentially 

/(r) — > -^exp(-0vr). (20) 

This means that the potential localizes the skyrmion exponentially. 

Of course we do not know whether such fields are the global minima of the energy. To check this 
we can perform a numerical simulation of (^) taking our derived field configurations ([ll| ) as initial 
conditions. We perturb them a little and then evolve them according to (^) with an extra dissipative 
term - jdt(j)a added to the right hand side of (||). We vary our perturbation and look at the 
configurations and their energies that our fields finally settle at. When the original configurations 
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are the global minima of the energy the perturbed fields evolved back to them; otherwise (for 
sufficiently large perturbations) the fields evolve to other static solutions and, hopefully, the global 
minima. Another possibility, in particular to check whether these solutions arc global minima, is 
to start with a more general field configuration corresponding to n skyrmions, say, n skyrmions 
located on a circle, and then evolve it using dissipative dynamics. In this case the circular set-up 
of n I-skyrmions with relative iso-orientation Sx = ^ was obtained by using the combination 
of 1-skyrmions determined by the hedghog field ( pi] ) and combined together as discussed in the 
appendix. 

The actual outcome {ie which fields are the global minima) depends on the form of the potential. 
In the next section we discuss the potentials that we have used in our simulations. 

4 Different potentials 

4.1 General comments 

A skyrmion is a S*^ ^ S'^ map and so as the field at spatial oo corresponds to (fi^ — +1 {ie the 
"North pole" of the field 5*^) and the skyrmion's position is at the point for which (p^ = —1. If 
V{4>3 = —I) does not vanish it costs energy for the field to take this value; hence in this case we 
expect the skyrmions to repel when they are brought "on top of each other" . Thus the minimal 
energy multiskyrmion configurations should be different when V{(j)3 = — 1) = and ^ 0. In 
fact, this was seen in the earlier studies; in the "new baby Skyrme" model {V{4>3) = 1 — 0|) the 
skyrmions lie on "top of each other" while for the "old baby Skyrme" model {V^cj)^) = 1 — (/is) they 
are separated from each other by finite distances. We have repeated these simulations and have 
extended them to a larger number of skyrmions. Below, we present our results: 

4.2 \/ = (1 + 03)4 

The first potential studied in the baby Skyrme model was V — [l + cji^Y Uli Hi)- This potential 
was chosen because the equations of motion have an analytic static solution (for n ~ \) oi the form 
W = X(x + iy), where A is related to 9s and By- The resultant soliton is only polynomially localised. 
Two such skyrmions repel each other and scatter at 90 degrees when sent towards each other with 
sufficiant speed. The model has further solutions corresponding to several solitons "on top of each 
other" ; but these solutions are unstable. When perturbed the solitons separate and move away to 
infinity. Thus this model have no stable multi-skyrmion solutions. 

4.3 V = l-(j)3 

This potential was studied in some detail in and |^ . As the value of the potential at r = {(j)^ = 
— 1) is V{(j)3) = 2 we expect some repulsion of the skyrmions when they are close together. Thus 
a field configuration of n skyrmions on top of each other is unstable. This was studied a little in 
1^ where it was shown that as a result of this repulsion the minimal energy multi-skyrmion field 
configurations produce nice lattice- like patterns. We have repeated these studies and extended them 
further to larger number of skyrmions. As there appear to be many solutions we have performed 
many simulations starting with different initial conditions and different perturbations. 

Our results show a rather complicated pattern of minimal energy field configurations. The 
results are presented in table 1. We use the coefficients of [|| ie add a factor ^ to the sigma model 
term, 63 = 0.25 and By = 0.1. 

All configurations are built out of 1, 2 and 3 skyrmions. In fig. 1 we show the energy densities 
of a pair and a triple. 
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Charge 


Energy 


Energy per skyrmion 


description 


1 


1.549 


1.549 


- 


2 


2.920 


1.460 


a pair 


3 


4.400 


1.466 


a triple 


4 


5.827 


1.457 


two pairs 


5 


7.321 


1.464 


triple + pair 


6 


8.781 


1.464 


6 individual 




8.731 


1.455 


three pairs 




8.790 


1.465 


two triples 


7 


10.248 


1.464 


pair+triple+pair 




10.230 


1.461 


triple+2 pairs 


8 


11.688 


1.461 


4 pairs 


9 


13.107 


1.456 


4 pairs + single 


10 


14.550 


1.455 


5 pairs 



Table 1: Multi-skyrmions of the old baby Skyrme model with V = 1 — (p^ 

A pair of skyrmions is very bound, more bound than a triple, and the pairs and triples are 
further bound for larger values of n. Thus the binding per skyrmion of a configuration of 4 and 5 
skyrmions is less that of a single pair or a single triple, respectively. 

The first more interesting case is of 6 skyrmions which seems to have several bound states. 
Of these, the state corresponding to 3 pairs is the most bound, the other two (2 triples and 6 
individual ones) are at best only local minima. For 7 skyrmions the lowest energy configuration 
corresponds to 2 pairs and a triple, which can be thought of a bound state of 4 and 3. The other 
state, corresponding to a triple sandwiched between two pairs has a higher energy and is at most 
a local minimum. In fig 2. we present energy densities of these two configurations. 

For higher n the situation becomes even more complicated; the lowest energy states always 
involve sets of pairs, but there seem to be other states involving triples and even singles. 

One question one can ask is whether the n = 2 state involves two skyrmions "on top of each 
other" or slightly displaced. When started from the original configuration of two skyrmions "on 
top of each other" slightly perturbed the system evolves a little, but this evolution is comparable 
to the original perturbation. Thus we believe the two skyrmions are slightly displaced, but this 
displacement is almost infinitesimal. 

4.4 y = 1 - 02 

The model with this potential was investigated in great detail by Weidig As both (j)^ — ±1 
correspond to the vacuum there is no energy argument which stops the skyrmions from "lying on 
top of each other". This is in fact what Weidig saw in his simulations; the states found by the 
shooting method (and so corresponding to many skyrmions "on top of each other" ) are the minima 
of the energy. Thus the energy densities of lowest energy multiskyrmion configurations have a 
ring-like structure. 
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Charge 


Energy 


Energy per skyrmion 


Break-up modes 


lonisation Energy 


1 


1.281 


1.281 


_ 


_ 


2 


2.332 


1.166 


1 + 1 


0.231 


3 


3.430 


1.143 


2 + 1 


0.184 


4 


4.540 


1.135 


2 + 2 


0.125 








3 + 1 


0.172 


5 


5.654 


1.131 


3 + 2 


0.107 








4 + 1 


0.166 


6 


6.771 


1.129 


3 + 3 


0.087 








4 + 2 


0.099 








5 + 1 


0.164 


7 


7.891 


1.127 


4 + 3 


0.078 








5 + 2 


0.095 








6 + 1 


0.162 


8 


9.011 


1.126 


4 + 4 


0.067 








5 + 3 


0.072 








6 + 2 


0.092 








7+1 


0.161 



Table 2: Multi-skyrmions of the baby Skyrme model with V = ^§(1 — (f>l) 



4.5 y = 01(1-02) 

The extra factor (^3 in this potential causes it to have an additional zero and so the vacuum structure 
is even richer. However, like for the 1 — potential, there are no energetic arguments which would 
disfavour the skyrmions lying "on top of each other". Indeed, this is what the simulations have 
shown; and we have checked this by starting with the hedgehog ansatz and n solitons on a circle. The 
simulations have shown that the hedghog ansatz configurations are the minimal-energy solutions. 
The shapes of energy density have radial symmetry and resemble the configurations of Weidig. The 
difference is in the energy levels (in this potential the energy per skyrmion decreases monotonically 
with n). Our results are presented in table 2. In this (and the next table) we have added a factor 
I to the sigma term and chose the 6 coefficients to be ^5 = 0.2 and 6v = 0.05. The values of the 
energies are divided by 47r. 

4.6 y=02(i_^3) 

Like in the previous case the extra factor (p'^ changes the nature of the solutions. In a way the effects 
due to this potential arc closer to those of the "old baby Skyrme model", in that Viip^ = —1) 7^ 0. 

We have found that for this potential only the hedghog ansatz with n = 1 has minimal-energy 
and that for 1 < n < 8 the skyrmions "on top of each other", when perturbed, evolve into 
configurations corresponding to n-skyrmions that lie on regular polygons. When starting from a 
circular set-up for this potential we have calculated the total energy of n 1-skyrmions as a function 
of the radius of the circle to get the initial state with minimum energy. This way the initial state 
settles down to the stable state sooner. Wc have confirmed this by performing also some simulations 
with different radi. We have also added to the initial configuration a non-symmetrical exponentially 
decaying perturbation to check whether the resultant multi-soliton states are stable, and whether. 
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Charge 


Energy 


Energy per skyrmion 


Break-up modes 


lonisation Energy 


1 


1.249 


1.249 


_ 


_ 


2 


2.372 


1.186 


1 + 1 


0.127 


3 


3.530 


1.177 


2 + 1 


0.091 


4 


4.693 


1.173 


2 + 2 


0.050 








3 + 1 


0.086 


5 


5.859 


1.172 


3 + 2 


0.043 








4 + 1 


0.084 


6 


6.926 


1.154 


3 + 3 


0.133 








4 + 2 


0.138 








5 + 1 


0.182 


7 


8.131 


1.162 


4 + 3 


0.092 








5 + 2 


0.100 








6 + 1 


0.045 


8 


9.300 


1.163 


4 + 4 


0.086 








5 + 3 


0.088 








6 + 2 


0.002 








7 + 1 


0.080 



Table 3: Multi-skyrmions of the baby Skyrme model with V = cf)'^{l — (p^) 



after dissipation, we reproduce the previous solutions. 

In fig. 3 we present picture of the energy density of the lowest energy field configurations 
involving 6 skyrmions. 

Our results show that this time (when we compare this case with V — 1 — ^3) the bindings are 
more comparable; in fact three skyrmions are more bound than two, and so the pattern is very 
different. It is interesting to note that the most bound system involves 6 skyrmions. 

Our results on the binding energies etc are presented in table 3. 



5 Comparison of potentials between 2 skyrmions 

We have also looked at the "potentials" between two skyrmions in our models. To do this we had 
to decide how to construct field configurations involving two skyrmions. We decided to do this as 
follows: We have taken f{r) for n — 1 and computed the field W from (^0|). Then we combined two 
such fields with Sx — and varied the distance between them computing the energy as a function 
of the distance. We examined the following ways of combining 2 skyrmions: 



W = Wi+W2 (21) 

W ^Wi~W2 (22) 

W ^Wi+W2-WiW2 (23) 

W = Wi-W2 + WiW2. (24) 

As could be expected, we have found that for very small r none of them gives reliable results. 
However, when two skyrmions are well separated the combination (|2^) gives us the lower energy 
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Charge 


Energy per skyrmion 


Energy per skyrmion 


1 


1.564 


1.564 


2 


1.405 


1.384 


3 


1.371 


1.344 


4 


1.358 


1.329 


5 


1.352 


1.322 


6 


1.349 


1.319 



Table 4: Comparison of energy per skyrmion between V = 1 



bl and V 



than ( pll) and so approximates the 2 skyrmion field more accurately. The figures presented below 
show our numerical results for the following potentials 



V^{1- 

y = (i- 



03 



(25) 
(26) 
(27) 
(28) 



As expected when the two skyrmions are far from each other the energy approaches twice the 
energy of one skyrmion. In our plots fig 5. we have indicated the energy of the hedgehog field 
with n — 2. Clearly, two skyrmions of this field are at r = but for better visualisation we have 
indicated it by a line. The low r part of each plot (with points indicated by A) overestimates the 
energy and cannot be trusted (at these points the skyrmions deform each other and are not given 
by ansaetze (21-24). However, as the real energies are lower we see that our results confirm again 
that for the potentials ( p6| ) and ( [27| ) the hedgehog fields are minimal-energy solutions but for (25) 
and ( p8| ) the minimal-energy solutions are different. The difference is more obvious in the (28) 
case as the plot indicates that the minimum of the energy is at r ~ 3.75. For (|25| ) it may seem 
that the hedgehog field with n = 2 is the minimal energy state but when we have evolved it the 
energy density expanded a bit thus showing that two baby skyrmions separate from each other a 
little. At the same time the energy has decreased and so we can conclude that for this potential 
the minimum of energy field configuration also differs from the hedgehog field. Incidentally the 9' s 
used to calculate the potentials (25) and ( p6| ) were chosen so that all four curves are compareable. 
Thus they are different from values used in table 1 and 4. 



6 Further Comments and Conclusions 

The extra coefficient 0| causes the potential to have a further minimum (at 03 =0). The question 
then arises of how this change effects the properties of the solutions. To perform such a comparison 
of multi-skyrmions of potentials (|2^) and (|2^) we have chosen the parameter 0s 's in both models 
to be such that the one skyrmion solutions of these models have the same energy. 

Table 4 shows the energy per skyrmion for the new baby skyrmion modelH and for the model 



with the potential (27) respectively. The results indicate that a further minimum in the potential 
increases the binding of the skyrmions {ie the multiskyrmions have lower energies). The same seems 
to be true for the other two potentials. 

We have also tried to understand the difference between the potentials (|28| ) and ( p5|) (in the 
(I25I ) the minimum of the potential is very close to r = while for ( ^ ) it is clearly much further 
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Figure 1: The total energy of 2 baby skyrmions as a function of r for different potentials (G, defined 
in (10) given by respectively, G(03) = 1, ^((/.a) = (1 + (^3), ^(.^3) = (1 + <P^)4>1 and G{4>^) = 4>l ) 
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out). A possible suggestion is that, in the latter case, the skyrmions move to a distance between 
them at which (j)^ — 0. Unfortunately this idea is not supported by the results of our numerical 
simulations. We have looked at the plots of the energy density and of the 03 at different distances 
for 2 baby skyrmions and have found (f)^ ^ 0. When we repeated this study for larger By (to have 
a more effective potential) the value of 03 decreased but was still nonzero. So the behaviour of 
skyrmions is more complicated and it depends on the properties of the potential in a more global 
way. 

We have also calculated the ionisation energies and various break-up modes of multiskyrmions 
for potentials (27) and (28). Comparing potentials (25) and (28) we see some further similarity (in 
each case the system of 6 skyrmions seems to be the most bound) and in each case, eg, the state of 
8 skyrmions needs very little energy to break up into (2+6). This is more pronounced in the case 
of the potential (28) for which the binding is stronger. 

In this paper we have looked at baby Skyrme models with more general potentials. We have 
found that as the skyrmions involve mappings between — > their properties depend crucially 
on whether the potential vanishes at the positions of the skyrmions (ie at 03 = —1) (we have 
assumed that at spatial infinity 03 = 1). When the potential vanishes the skyrmions "lie on top of 
each other", when it does not - they separate and form interesting lattice like patterns. The shape 
of these patterns depends on the details of the potential. The same holds for the binding energies 
of skyrmions in all models. When the potential is more complicated (ie it has further zero) the 
skyrmions are more bound and in their patterns are more spread out - however, the actual positions 
and distortions depend on the details of the potential. 

When the skyrmions are spread out the system of skyrmions has many local minima with some, 
larger or smaller, potential barriers between them. Thus for instance, for = 1 — 03, a system of 
6 skyrmions has at least 3 local minima, and depending on the initial configuration the system can 
land in any of them. This is not unexpected and it suggests that similar models in, physically more 
relevant, 3 spatial dimensions may also have many local minima. Thus, the problem of finding 
multiskyrmion solutions of models in higher dimensions is clearly very complicated. 
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7 Appendix 

Here we make some remarks about our numerical procedures 

Profile functions We have used the shooting method to determine /(r) and have integrated 
( [Tsl ) by a fourth-order Runge-Kutta method for any n. To avoid a singularity at r = we have 
considered r ~ e with e small. We have used the formulae (17) and ( p^ ) with 10000 lattice points 
and with the spacing dr — 0.003. 
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Hedgehog static solutions In most of our simulations we have used a 201 x 201 lattice with 
lattice spacing Sx = dy = 0.3. For each (xi,yj) f{xi,yj) was determined by a linear interpolation 
using the values determined by the shooting method. Given / we calculated cj) using the hedgehog 
field expression (nih 



A linear superposition for static solutions with n > 1 In an alternative formulation 
we have used a single complex field W, which is related to (f> by 

W + W* , W~W* , 1-WW* 
4>i = 02 = I 03 = (29) 

where * denotes the complex conjugation. Hence the complex field W when expressed in terms of 
the profile function /(r) and d takes the form 

W = tan f ^) e-'"^ (30) 



Static initial field configurations with Q = n were formed by a linear superposition. When the 
baby skyrmions are far from each other (but not too close to the boundary) we can construct a 
configuration with charge n from W with n = 1 by a linear superposition 

W{x,y) ^^Waix- xa,y- ya), (31) 

a 

where {xa,ya) is the location of the centre of the ath skyrmion and Wa is the field ( |30| ) with n = 1 
and appropriate f{r). We have used this method of constructing initial configurations in some of 
our simulations. Namely, we have placed n baby skyrmions at equal distances from the origin with 
relative phase shifts Sx = — between them (for maximal attraction) and then used (^9|) to get $. 



Time evolution of the static solutions We integrated the equation of motion for each 
component of cj) independently. In this manner, during the simulations the field would gradually 
move away from the unit sphere S^. To correct this, every few iterations, we kept on rescaling the 
field as follows: 

" = (32) 



and 

9t0a > dt<j)a - t (t>a- (33) 

Cf) ■ (j) 

Another problem we have had to face involved using a finite lattice. Thus we had to make sure 
that the boundary effects did not alter our results. To be certain of this we varied lattice spacing 
and the number of lattice points. 

When the fields were not at the minimum of the energy we allowed them to flow to this minimum 
- reducing the energy by a damping term. 

dtt4>a = K^bTb {$, dt$, di^'j - ^dt(j)a, (34) 
where 7 is the damping coefficient. We set 7 to 0.1. This term takes the energy out of system. 
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We have used a fourth-order Runge-Kutta method to simulate the evolution of the field in time 
working in double precision. The time step was 6t = 0.1 throughout. We have performed our 
simulation using the (j) formulation with the derivatives replaced by finite difference^ as explained 
in Q, see also [||). 
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